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1.0 INTRODUCTION 


This report provides a description of two methods and their computer 
programs for calculating viscous and inviscid flows past axisymmetric bodies. 
The first is a boundary-layer method for both laminar and turbulent flows. 
The procedure embodies a finite-difference representation of the boundary- 
layer equations and uses an eddy-viscosity model to represent turbulent 
flows. The details of the eddy-viscosity model and an indication of its 
range of applicability for flows past surfaces with heat and mass transfer is 
provided by Cebeci and Smith u ._ The two-point finite-difference method is 
described by Cebeci and Bradshaw 2) Section 2 describes the method and states 
the boundary-layer equations, corresponding boundary conditions, fluid 
properties, turbulence model, solution procedure and the boundary-layer 
parameters. Section 3 describes the computer program and states the input 
and output of the computer program together with three sample calculations 
to demonstrate the capability and use of the computer program. 


The second method is an inviscid-flow procedure developed by James 2). 
Solutions to the potential equation are obtained by an orthogonal expansion 
procedure. For the cases to be presented, the accuracy is equivalent to 
higher-order Neumann solutions, but at a considerable savings of computer 
time. When used in the direct mode, this computer code yields the pressure 
distribution (edge velocities for boundary-layer calculations) for any 
specified body shape. An additional feature of the inviscid procedure is 
its capability to be used in an inverse (or design) mode. Here, the desired 
velocity on the body is specified, and the procedure will yield the body 
shape consistent with this velocity. Details of the development of the 
method are given in Section 4. Section 5 describes the use of the computer 
program, and demonstrates its capabilities on five different body shapes. 


2.0 BOUNDARY-LAYER METHOD 


2.1 Boundary-Layer Equations 


The governing boundary-layer equations for steady axisymmetric 
compressible laminar and turbulent boundary layers are the continuity, 
momentum, and energy equations. These equations and their boundary 
conditions for the coordinate system shown in figure 1 are: 


Continuity 
3 (our) + = (svr) = 0 (2.1) 
ax ‘? ay *° 
Momentum 
oe Mey, Geald [of a parr] (2.2) 
) _  e y e dx r oy wae : : 
Energy 
ays ake ta fly at pry (fy —L) yu 
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Figure 1. The notation and the coordinate system for axisymmetric flows. 
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ay Ny 
y>e  uru(x), H+ H, (2.4b) 
Here 
pV = pV + pv" 
r=r, + y cos¢ (2.5a) 
dr 
6 = tan”! 2 (2.5b) 


We note that with r being a function of x and y, these equations 
account for the transverse curvature effect which becom:s important when 
the boundary-layer thickness is of the same order of magnitude as the body 


radius ro: 


The solution of the system given by (2.1) - (2.4) requires fluid 
properties and closure assumptions for the Reynolds stresses, -pu'v’ and 
-ov'H'. In this report we consider both air and water and discuss their 
fluid properties in Section 2.2. To satisfy the closure assumptions, we 
use the eddy-diffusivity concept and represent them by the algebraic eddy 
viscosity formulation of Cebeci and Smith‘ °. The relevant equations for 
the eddy-viscosity formulas are discussed in Section 2.3. 


2.2 Fluid Properties 


The solution of the boundary-layer equations recuire the specification 
of fluid properties. In our study we consider air, pure water and sea 
water as the working fluid and define their properties below. Air is treated 
as a perfect gas, and the fluid properties u and p are assumed to be 
functions of temperature only. The specific heat of air at constant pres- 
sure, namely, Cys is assumed to be constant and equal to 6006 ft2/secR. 
The viscosity is obtained from Sutherland's law expressed as 


2.270 10° ; “es ta i (2.6) 
we. X ° 
8 Ts 108.6 — a 


The density-temperature relation is obtained from the equation of state 
and from the assumption that the static pressure remains constant within 
the boundary layer, 


1b -sec? 
= er — (2.7) 
o* WaT ae 


The Prandtl] number of air is an input to the computer program and is assumed 
to be constant. 


The fluid properties of pure water are given by: 


u = 3.74614 x 107° / (35.155539 — 106.9718715T,, + 107.7720376T° 


— 40.5953074T3 + 5.6391948T?) (2.8) 
p = 1.94 (0.803928 + 0.4615901T,, — 0.2869774T% + 0.0234689T2) (2.9) 
cy = 2522.98 (1.483689 — 0.8072501T, + 0.3289607T") (2.10) 
Pr = 13.66 / (73.376906 — 208.7474538T,, + 197.7604676T~ — 68.8626188T° 
+ 7.4779458Ts) (2.11) 
k = 6.917333 x 10 (-1.9410583 + 5.2220185T, ~ 2.693322T? 
+ 0.4176176T2) (2.12) 
The fluid properties of sea water are given by: 
u = 3.9258 x 107° / (-289.99547 + 1101.820826T, — 1577.2442T2 
+ 1003.8069T? — 237.3881T)) (2.13) 
p = 1.9947 (0.66285108 + 0.69246952T,, — 0.3553206T*) (2.14) 
c, = 23859.32 (1.774152 — 1.518056T, + 0.7439018T°) (2.15) 
Pr = 13.565 / (-686.60487 + 2653.55416T,, — 3851.34007TS + 2482.21458T° 
— 596.8245T") (2.16) 
k = 6.90892 x 10° (12.28011 — 35.655012T, + 36.562345T? 


Yr 
fe 12.187448T}) (2.17) 


Here Li is a temprature ratio defined by 


(2.18) 


The units of u,o,¢,,k and T in (2.8) - (2.17) are Ib,-sec/ft”, 
slugs/ft?, ft2/secR, (1b,-ft)/sec ft°R, °R, respectively. 


2.3 Turbulence Model 


To satisfy the closure conditions we use the eddy viscosity (e,) and 
eddy-conductivity (e,) concepts and define them by 


ou sTor — oH 
PED, ay? -ov'H' = pen ay (2.19) 


and relate En and e, toa turbulent Prandtl] number Pry by 


(2.20) 


According to the eddy-viscosity formulation of Cebeci and smith! ), 
the turbulent boundary layer is divided into two regions, called inner and 
outer regions, and eddy-viscosity formulas are defined by separate formulas 
in each region. They are: 


_ 2 [au 
=U? [tlre (egy < (eq)y 


(2.21a) 


= 0.0168 J — u)dy Vee (endo > ted 


F (2.21b) 


L = 0.4y[1 — exp(-y/A)] (2.22a) 


A = 26v/N ur! (o/o,)'?, u, = (r/o) !/? (2.22b) 


Tt 


For flows with no mass transfer the parameter N_ is 
1/2 
We\ (Pe \2 
N= |] —11.8 i p 
He /\Pw 
When there is mass transfer 


0.\2 _+ u 1/2 
my 8 ae a ‘) ae 
N | (; ) ; [ exp (11.8 7 + exp (1.8 : ‘I 
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where 


Vv 
p = BF ne i vee ot (2.22e) 


The parameter Tee allows for the finite length of the transitional region 
between a laminar and a turbulent flow and is given by 


x Xx 
ee by dx dx 
Yep 2 1 — exp |-6(r,) I a 4 - (2.224) 
Xer Xtr 


Here G is a spot-formation-rate parameter 


3 
fafa late Ae 
Ye 5 tr 


According to the studies conducted by Cebeci and Smith, a constant value of 


turbulent Prandtl number is quite satisfactory for the type of flows to be 
considered in this study. Therefore, we assume Pry = 0.90. 


2.4 Mangler-Falkner-Skan Transformation 

Before we solve the boundary-layer equations for both laminar and tur- 
bulent flows, we uSe a combination of Mangler and Falkner-Skan transforma- 
tion. For this purpose we define new independent variables xX and n_ by 


« Pi 2 
dx = (2) dx (2.23a) 
u 1/2 
Bont ean r 
dn = on o(E)ay (2.23b) 
and by introducing a stream function that satisfies the continuity equation 
(2.1) 
—“ (2.23 
- - .23c) 
wv (o,ugUeX) Lf(x,n) 
a so that 
KI onus 
Hy pur = = , ovr = - + a (2.24) 
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With this transformation it can be shown that the momentum and energy equa- 
tions and their boundary conditions can be written as 


Momentum 


ax ax 


ag : PLM A YS 


(b,f") + m ff" + mc — (#*)?] matt = % (r af" _ ¢" *t)(2.25) 


Energy 
(baE") + (byft#") + (mf —m,)E! = x (r _ =) (2.26) 
aX ox 
n=0 f=f'=0 ag, + o£! = given (2.27a) 
n>n f' +1 E+1 (2.27b) 


Here primes denote differentiation with respect to n and 


= 2 + = py = Pe ~ = em 
b, cil +t)(1 +e), Cc fale” oar es 
2 (2.28) 
ml hse > es 2 EE at 
bo =e) Pr) i b3 = Pr (1 + t) ( ™ Pr, “) 
- - du 
ee een. oe 
Ry Sh gt Se (ogi) we eae 
@ ax 
0, Vv u x 
my = ww pil? (=) : Rs = as (2.29) 
Pole 0 Ve 
2(%/L) cosél, 1/2 n 
B® <P t ee ee > 12 fican 
(r/L)" (Rs) 
H 
fi=u, Eso (2.30) 
Ue He 


Similarly the eddy-viscosity formulas given in Section 2.3 can be expressed 
in transformed variables and can be written as 


(en), = 0.16 c@ (2) (re) (ty pel? (+ ¢)e? | — exp (. £.- fel) If" [Yen 


cos $ 
(2.31a) 


Poll —F 
c(l -f' 
f aye 


+) Heh (t \ 1 
(c*), = 0.0168C (-2) (t) Ru/ tr (2.31b) 
0 
Here 
r r_\3/2 
0 = N_ 3/4 0 L -3 2 “ 
PREY ge aw 
1/2 
a 4 Mw) ("0\ 9-1/2 en 
fe aig se usu, (“*)(c2) 5 (Ft) (2.32) 
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2.5 Solution Procedure 


We use the numerical method described in ref.4 to solve the governing 
equations presented in Section 2.1 and 2.2. This is an efficient two-point 
finite-difference method developed by Keller and cebeci (4) and extensively 
used by Cebeci for two-dimensional and three-dimensional flows (see, for 
example, ref. 2). A detailed description is presented in ref. 2 and is not 
repeated here. 


One of the advantages of the present numerical method is that nonuniform 
net spacings can be used in the x-direction as well as across the boundary 
layer. In the latter case, the nonuniform grid is a geometric progression 
with the property that the ratio of lengths of any two adjacent intervals is 
a constant; that is, An, = Kans_}- The distance to the j-th line is given 
by the following formula: 


any = an, (« wie =1) KI (2.33) 


There are two parameters in (2.33): Any» the length of the first step, and 
K, the ratio of two successive steps. The total number of points J can be 
calculated from the following formula: 


In[T + (K —1) (n/n, )] 
Ce (2.34) 


In the computer program which embodies the present solution method, Any and 
K are specified. For moderate Reynolds numbers, say Ry up to 107, typical 
values of Any and K are 0.01 and 1.14, respectively. In general, approx- 
imately 50 grid nodes across the boundary layer are sufficient to represent 

laminar and turbulent boundary layers. On the other hand, at higher Reynolds 


Yin et WAGs bh Be SORTER 


ne ne 


Te 


numbers it is desirable and necessary to choose Any equal to 0.001 and K 
to 1.26, respectively. 


In the present method, the flow is laminar at the stagnation point and 
can become turbulent at any specified station greater than second x-station 
(NX > 2). The present boundary-layer program does not compute the transi- 
tion location and requires it to be input. 


2.6 Boundary-Layer Parameters 


The output of the computer program includes profiles of stream function, 
velocity and shear stress, i.e., f, f' and f" as a function of the sim- 
ilarity variable n and the physical distance y at calculation station x. 
Profiles for the total enthalpy E distribution and its derivative E' are 
also presented. The calculated boundary-layer parameters such as the local 
skin-friction coefficient Ces the total drag, the various boundary-layer 
thicknesses and shape factors are summarized at the end of the output data. 
Expressed in either physical or transformed coordinates, these parameters are: 


The local skin-friction coefficient: 


at, Cy (To 
oe io) ¢e (2.35) 
ee X 


The drag coefficient based on total frictional resistance D is: 


X <— 2 
—. 2D... 2 29 e\/-e 
a aah f cet (:2\() dx (2.36) 


Here L is any convenient reference length defined by the user and x is 
the axial distance from the nose of the body. The subscript t refers to 
the aft end of the body. To calculate the total drag of the body, use is 
made of Young's formula, which can be written for compressible flow with 
average values in the wake as 


2 
2(6n)+ ( 2.5 uphetS*0.44 0/2 
Cy = eee 2) (2) (2.37) 


t t 


oq = 2n i r (EV(i-V( - ile = = fr — f')dn (2.38) 


RS 


and h is the shape factor based on the ratio of the displacement area oh 


to the momentum area 6y3 SK being defined by 


aoe fe (1. - mr a} y= st (c£ — +1) dn (2.39) 


Sometimes it is useful to examine integral thicknesses based on the profiles 
alone; that is, thin two-dimensional equivalents: 


of S( “a ir (FE Alisa dy (2.40) 


par 


p tiie Bike). ee FE 
of Ae ( ah + tre (i) Sf Tha & aie 
0 x 0 


H.. +. (2.42) 


The height of the displacement surface, s§ for axisymmetric flows 
with transverse curvature term included is defined by: 


2 X r 1 - 
Zehr, + (8§)° cose = 2 f r ( - <A ay = eset Jan = — (2.43) 
ee X e 

0 0 
or 

* r. 4 cosé V2 

Soap | ee eu 

ro 

Note that for "> 83 the above expression gives 83 > 85 as expected. 


The boundary-layer thickness 6 jis obtained by interpolation for y 
where u/u, = 0.995. 


3.0 DESCRIPTION OF THE BOUNDARY-LAYER COMPUTER PROGRAM 


3.1 Input 

Given an axisymmetric body at zero: incidence and the associated 
inviscid velocity distribution on the surface of the body, this program 
calculates the development of laminar and turbulent boundary layers with air 
or water as working medium. Surface boundary conditions allow for arbitrary 
wall mass transfer, wall temperature or heat transfer rate. 


For air the external inviscid flow is treated as compressible isentropic 
flow. The external temperature and pressure distributions are calculated 
from the following formulas: 


T u 2 

6 A ee Be =) (3.1) 
a 2 o |\u, 

ax y/y-1 

see -re tel) ai (3.2) 


Here y is the ratio of the specific heats at constant pressure and volume; 
its value is set equal to 1.4. The total enthalpy at the boundary-layer edge 
is calculated from 


2 
Ue 
H, = ¢.T, + 5 (3.3) 


and density is calculated from the perfect gas formula (2.7). The kinematic 
viscosity Up is calculated from (2.6). 


For flows with water, either pure water or sea water, the fluid is 
assumed to be incompressible. That is, there is no Mach number effect but 
all fluid properties are allowed to be functions of temperature. In addi- 
tion, the specific heat at constant pressure is assumed constant and equal 
to the average of the wall and freestream values. This is a reasonable 
assumption for the temperature levels and differences occurring in practice. 
As a result, the calculation of temperature from the total enthalpy is 
greatly simplified. Fluid properties for water are calculated by (2.8) to 
(2.17). 


11 


Essentially the input to the computer program consists of five types 
of cards. Card 1, shown below, contains the title of the flow problem 
under consideration. 


otoaben 
Weae 


nL af3]<]slolz]e] heheh shehep fee shed pape 


TITLE 


bebo bepabaorbebord 


ddccce 


a Ob 8g ee 8 Pk 
Load Sheet for Card 1 
The input is punched as 80-column alphanumeric field. 


Card 2 requires the following information to be specified. The input 
is in 613 format. 


HEBAARBROERIELERIE 
xt | wre] kN |oSFO| nxy [IWALL 
ava Tas he OE 


Load Sheet for Card 2 


NXT total number of x-stations to be calculated NXT>3 

NTR station number at which turbulent flow calculations start. 
NTR>3 must be observed. 

KN this flag controls the type of fluid to be used in calculations 
= Air 
= 2 Pure Water 
= 3 Sea Water 

ISFD this flag controls the surface distance calculations. Usually 


we input pairs of the body coordinates at the x-stations where 

the boundary layer is to be calculated. But sometimes it is 

advantageous to use a different set of body coordinate data or 

input surface distance directly. 

O surface distance calculated directly at x-stations. 

1 surface distance calculated from the separate body coordinate 
input and interpolated to desired x-stations. 

2 surface distance input at x-stations. 


number of body coordinate pairs for separate surface distance 
calculations. Used only if ISFD=1. 


flag that controls the wall boundary conditions for the energy 
equation. 

1 wall temperature specified, °R. 2 

2 wall heat transfer specified, Btu/hr/ft 


Card 3 contains the following information to be specified. The input is in 
8F10.0 format. 


AeT eT sTelelol oleh ahhed sel de eee ded ebeheh hsp preheat feehegeehopehpe 
SSR Ce ARN RTM 
bbs oh duhadoraberer Po Poy Peyees 


a ee ae 


Load Sheet for Card 3 


transformed boundary-layer thickness, n_. This is for the first 
station. A value of 6 is usually sufficient. For NX>1, the 
boundary-layer thickness is computed internally. 

DETA(1) initial An-spacing at the wall, Any. 

VGP variable grid parameter, K. 

UFRS freestream velocity u_, (ft/sec) 


TFRS freestream static temperature (°R) 


RE unit Reynolds number based on freestream conditions, Re /ft = u./v_. 


CL reference length L. Converts nondimensional body coordinates 
into dimensional form, expressed in feet. 


PR Prandtl] number for air. Can be ignored for calculations with 
water. 


The fourth and the following cards contain the information for surface 
distance calculation, the external velocity, surface temperature or heat 
transfer, and mass transfer at the wall. There will be NXT cards of this 
type. A sample load sheet for this case is shown below. The format is 
5F10.0. 


fala] sfel [el hota heheh ohhh ef add ade bop pepsi 
DE Broker Sirs Bors ote lh 


Load Sheet for Card 5 


If ISFD=Q (surface distance calculated directly) 

XC (1) nondimensional axial distance (X/L) from the nose of the body. 
RC(1I) nondimensional body radius (ro/t) 

If ISFD=1 (surface distance interpolated) 

XC(1) nondimensional axial distance (X/L) for desired x-stations. 
RC(I) can be left blank. 

If ISFD=2 (surface distance input) 


XC(I) nondimensional surface distance (x/L) input. 

RC(1) nondimensional body radius (r_/L) input at stations corres- 
ponding to (x/L) input. ’ 

UE(I) nomdimensional velocity distribution (u,/u,) input at x-stations. 

HTP(I) dimensional wal] temperature or heat transfer. 


If IWALL=1 the wall temperature (Ty) is input, °R. 
If IWALL=2 the heat transfer rate at the wall (qy) is input; 
Btu/hr/ft2. 


Bv(I) nondimensional mass-transfer ratio at the wall, pwYy/oel,- For 
no mass-transfer, leave blank. 


If ISFD=1 (surface distance interpolated) additional data is needed 
for the separate body definition cards. These are in the form of axial 
distance (X/L) and body radius (r,/t) pairs. There are NXY cards of 
this type. A sample load sheet for this case is shown below. The format 


is 2F10.0. 


HABQAABOOEELELERLEES 
ee ele 
Hh ph Perens 


Load Sheet for Body Coordinates 


XS(1) nondimensional axial distance (X/L) from body nose. 
RS(I) nondimensional body radius (ro/t). 
3.2 Output 


The output of the computer program includes printout of the data used 
in boundary-layer calculations as well as tables of profile data at all 
x-stations. At the end of the printout is a summary of the calculated 
boundary-layer parameter. Boundary-layer calculation data is presented 
under the headings of CASE DATA and STATION DATA: 


CASE DATA: Data is printed out in two groups. The first group contains 
the flags and control variables specified in the input. The second group per- 
tains to the input and calculated freestream data. 


UFRS freestream velocity (u,), ft/sec. 

TFRS freestream static temperature (T_), °R. 

PFRS freestream static pressure (P_), 1b/ft2. 

RHOFRS freestream density (p,), slugs/ft?. 

PRFRS Prandtl] number for air. 

RMUI freestream dynamic viscosity (u)> Ib-sec/ft® 

RCI unit Reynolds number based on freestream (u_/v.)> rt! 


Note that for air there is additional printout for the freestream Mach number 
(M_) and the freestream total temperature (T7), R. 


STATION DATA: refers to a table with the following columns. 
NX station number of the x-station 
X calculated or input surface distance x, ft 
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The profile data is presented in tables identified by a printout line con- 
taining the station number NX, the surface distance x, and the surface 

The next few lines before the table proper 
The columns in the table contain 


distance Mangler variable XB(x). 
refer to the iterations during calculations. 


input or interpolated nondimensional body radius (ro/L) 
inviscid velocity at an x-station, ft/sec. 

input heat-transfer data 

IWALL=1, wall temperature (Ty) distribution, °R. 


=2, wall heat transfer (qy) distribution, Btu/hr/ft* 


input mass transfer (po. Vv. /p_u_) distribution 
ww ee 


calculated pressure gradient parameter (m, = x/U, du,/dx) 


the following: 


ETA 


TEMP 


OUTPUT SUMMARY : 


point number: profiles are plotted from the wal] outward. 


all points are printed. 
nondimensional boundary-layer variable n. 


nondimensional velocity in the boundary layer, f' = u/u,. 


derivative of U with respect to n, f". 
nondimensional total enthalpy ratio, E = H/H, « 
derivative of W with respect to n, E’. 
static temperature in the boundary layer, °R. 
normal distance from the surface y, ft 


with two rows per NX-station. 


NX 


First row 


the x-station number 

surface distance (x), ft 

inviscid velocity (u,)> ft/sec 

boundary-layer thickness 6, y-value where u/u, = 0.995. 
momentum area 6,5 rt? 

two-dimensional momentum thickness 85» ft 

local skin-friction coefficient, Ce = 2, Joga 


This table at the very end of the output contains 
some input data and the calculated boundary-layer parameters in 9 columns 


Not 


CD 


TW 


Second 


X/C 


CT 
DISPLT 
QW 


total drag coefficient of the body Cp, calculated by Young's 
formula. Reference areais the input reference length (L) 
squared. 

wall temperature (T,)> R. 


Row 

the axial nondimensional distance from the nose of the body. If 
surface distance is input (ISFD=2) the input nondimensional 

x is printed out. 

local Reynolds number Re = UAX/ Ve 

displacement area Sh Ft? 

shape factor used in total drag formula Hy = Sy/6y 
two-dimensional shape factor H, = 53/85 

the average integrated skin-friction coefficient of the body, Ce 
displacement thickness for axisymmetric body 83, 2 

heat transfer rate on the body surface i “* aT/ay, Btu/hr/ft? 


3.3 Sample Calculations 


The present computer program can be used to determine the properties of 
axisymmetric laminar and turbulent boundary layers in air or water over smooth 
surfaces with wall mass transfer, arbitrary wall temperature or heat transfer. 
Flows in air are treated as compressible and flows in water as incompressible, 
but with allowance for temperature-dependent fluid properties. In order to 
demonstrate its capability and input requirements, we present these sample 
calculations for the same body representing flows in different fluids and at 
different wall boundary conditions. 


The body called 5A has a fineness ratio of 6 2/3 and its profile is 
formed by a modified NACA 4-digit thickness distribution with maximum thick- 
ness at 45 percent. Its nose radius and trailing-edge angle are standard. 

The body length is 20 ft and all runs were made at a length Reynolds number 
of 40 x 10°, This combined with a freestream temperature of 520°R corres- 
ponds to a freestream velocity of 570 fps in air and 35 fps in water. Transi- 
tion is fixed at about 21 percent of chord (station 15) for all cases. The 
external inviscid flow was calculated by the Douglas Neumann program. It 
should be noted that no compressibility corrections were applied to the calcu- 
lation in air although the freestream Mach number was 0.5. 


The sample case for flow in air was computed with adiabatic wall boundary 
conditions; that is, a. * 0, and suction was applied to the rear portion of 
the body from approximately 56 percent chord-point onward. The input data, 
typical turbulent boundary-layer profiles with suction at 81-percent chord 
(station 39), and the summary of calculated data are shown below. It should 
be noted that in this, as well as in two other cases, boundary layers sepa- 
rate at the last calculation station. This is expected since the inviscid 
flow at the aft end of the body approaches a stagnation point: If the pres- 
sure changes due to viscous interaction are not accounted for, continued 
boundary-layer calculations will result in inevitable separation. 
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The second sample case is for flow in pure water (KN=2) where the nose 
portion of the body under the laminar boundary layer was heated. A moderate 
heat flux amounting to 1 Btu/sec/ft2 was specified over most of the area, 
with less near the nose where a strong favorable pressure gradient exists. 
The input and output data for this second sample case are shown below. 
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The third sample case is similar to the second one except that the 
working fluid is sea water (KN=3) and calculations for the heated portion 
of the body were done with the wall temperature specified at a constant 
540°R. The surface temperature was assumed to vary linearly from its speci- 
fied value on the heated portion to the freestream value on the rest of the 
body. Summary of the input data is presented below. Profile data at sta- | 
tion 14 near 19-percent chord location is added to show the effect of heating 
on a laminar profile. Although the local pressure gradient is slightly | 
adverse, the two-dimensional shape parameter Ho from the summary page shows 
that the profile effectively corresponds to one obtainable in a slightly 


favorable pressure gradient with no heating. 
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4. STANDARD AND INVERSE POTENTIAL-FLOW METHOD 


The standard and inverse potential-flow solution procedure is based on 
the analytical method of James '3), This method, originally developed for 
shapes having either a cusped or blunt trailing edge, has been extended to 
include trailing edges of arbitrary angle. It is based on a mapping of the 
body shape to a unit circle. This mapping can be performed either analytically 
if the mapping function is known, or numerically through a program which 
determines the Fourier coefficients of the mapping. Once this is accomplished 
the surface coordinate of the body is expressed as an expansion in the circle 
plane, which allows the solution to the governing equation to be expressed as 
sequence of solutions to recursively generated equations. This permits an 
arbitrary number of terms to be generated depending on the accuracy desired. 
The details of this procedure are outlined below, and a description of the 
computer program is included in Section 5.0. 


4.1 Equations of Motion 

This section deals with the development of a general analysis method 
for axisymmetric flow very briefly. For more details, the reader should 
consult ref.. 3. 


For an axisymmetric body lying in the (x,y)-plane such that the free- 
stream at infinity is parallel to the axis of symmetry (x), the velocity 
potential (¢) satisfies the well-known differential equation 


sled) By 63) 9 o 


with the additional conditions that » ~x in the far field and that the 
derivative of 4 normal to the profile vanishes on the profile. It is 
quite standard to seek orthogonal curvilinear transformations x = X(EysEo)» 
i y = y(&y E>) in order to either simplify (4.1) or the profile representa- 
tion (thus making the normal derivative condition more tractable). Such a 
transformation leads to 


: h h 

F) 2. 36 p) 1 | a¢ 

3c, \he * bec) * SEL lh > ; 

| os (i : = ‘ “re ( 2 ; t-) s sans 


where hy and ho are the usual metrics, viz. 
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2 2 oe 
hy = + (24), hy (2) + (2) 
] ra d6y 2 Eo dEo 
However, if the transformation is in particular a conformal mapping from the 


ey + 1€5-plane to the x + iy plane, then the Cauchy-Riemann equations 
hold, that is 


and therefore, hy = h, which reduces (4.2) to 


a we 


Even though (4.3) is formally no simpler than (4.1), it is possible to find 
general mappings which greatly simplify the surface boundary condition by 
relating the given profile to one of a number of canonical forms in an 
auxiliary (€) &>)-plane. Specifically, the class of mappings which map 
airfoil-like profiles into the unit circle have been found to be of great 
utility in two dimensions and so are considered further for this application. 


4.2 Unit Circle Mapping 


According to Riemann's fundamental mapping theorem, any profile whose 
contour encloses a singly-connected domain (in the usual engineering sense 
of, e.g., Woods 6 ) can be mapped into the unit circle. However, for 
airfoils and axisymmetric bodies of interest to engineering, the types of 
mapping of sufficient generality and simplicity to be useful are quite 
constrained. The role of these constraints has been discussed by James 
and need not be elaborated in detail here. Briefly, the desired character- 
istics are that: (a) the mapping should reflect the existence of a finite 
trailing-edge angle discontinuity (7 —t), (b) the mapping derivative must 
+1 as |z| +, (c) the profile must be closed, and (d) the mapping de- 
rivative must have neither zeros nor singularities on* or outside the unit 
circle. 


(7) 


*In the general theory of mapping methods, certain boundary singularities may 
4 arise, but they are not of interest in this application. 
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Let ¢=€& + in be the unit circle plane then it is easy to show that 
a mapping which conforms to (a) - (d) is of the type 


- (0-7) ao ee 


where g(c) is an analytic function with neither zeros nor singularities 
outside C (the unit circle) and which +1 as |z| +. Many different 
forms of g(z) can be used to generate profiles, but for the present 
purpose the more general representation of g(¢) as an expansion 


a 
a(c) = 1+ av, By (4.5) 


. z 


Note that the first coefficient ay is given by a, = 1 —t/m because of 
the closure condition. 


Now y is required in (4.3) and is described in (4.4) as the imaginary 
part of z. It can be shown that z is governed by a factor of the form 
c(1 =tey ", and direct inclusion of this into the differential equa- 
tion would be desirable. However, an analytical solution under these cir- 
cumstances has not been found. An alternative simpler procedure is to 
expand (4.4) and (4.5) to give 


se SA RR Fae (4.6) 


for all values of +. The convergence of (4.6) is guaranteed when +t = 17 

or t= 0 because of the convergence of the a, coefficients — which is 

in turn guaranteed by the theory (see ref. 7). For arbitrary + (0< 1 < 7) 
it might be expected that poor convergence could be a feature very close to 
the trailing edge, but experience has shown (see ref. 3 and the results 
quoted later) that this is not an important consideration with the kinds 

of bodies and term numbers of interest to this work. Moreover, this is 

not a subject to be pursued in detail in this report. Accepting that (4.6) 
is sufficient, it is clear that the B, coefficients can be calculated 
readily from the ay, coefficients and the expansion for the singular factor 
in (4.4). On the other hand, if only a body shape is given, there are mapping 
programs available to find the ay coefficients. The programs associated with 
this theory contain a subroutine option of this type. 
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With regard to the choice of variables in mapping from the unit circle 
it should be noted that the natural group ¢ = re’ is such that z is not 
an analytic function of r+ iw so that it is necessary to work with ¢ = e°, 
o=Atin, t= e*. Then for axisymmetric bodies, the reflection symmetry 


ensures that all the ans B, coefficients are real giving, from (4.6), 


y =e sinw + Be sinw + Be“ sindw +... (4.7) 
to be used in the differential equation 
ox & HH) + 5 HH) - (4.8) 
with the conditions 
¢% e* cosw when A+ oe (4.9) 
89. 9 when d= 0 (4.10) 


Although this is a properly posed problem which has so far defied analytic 
solution, the reduction presented here is by no means exhaustive; it is 
merely the simplest. 


4.3 Method of Solution 


The structure of y (eq. 4.8) suggests a scheme of solution in which 
¢ iS computed from a sequence of potentials $n sensitive to sin nw at 
each stage. An endless variety of such schemes exist, @.g., Kaplan | 8) 
His procedure is both elegant and complicated. Too complicated, in fact, 
to proceed either algebraically or automatically beyond the sin3w term. 
From this,one rather obvious inference is that if the scheme is to be 
reducible to automatic computation of coefficients (as required to substanti- 
ate any claim of generality or analyticity) then each successive potential 
or function must not be governed by another partial-differential equation — 
and not even by an ordinary differential equation unless it can be solved 
directly in terms of known functions. With this understanding the natural 
choice for » (an even periodic function) would be 


> = F(a) + £,(A) cosw + f(r) cos2w +... 
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or perhaps, using & = cOSw, = e*, 


6 = gy(n) + g)(nde + gp(nde? +... 


However, these lead to systems of equations with strong coupling between 
the n-th function and higher order functions not yet known. 


It turns out that the format which works in the sense that all coupling 
between equations is to lower order known functions is 


@ = nFo(e) + LF (e) + wz Fale) + oe. (4.11) 


and this is the central feature of the whole theory. To go with (4.11) the 
sinw factor in y is removed reexpressing (4.7) as y = sinw - Y(E,n) 
where 


B, 8, 
Y(é,n) = n+ —* ls a 


and U, Ce) is the Chebyshev polynomial of second kind. At the same time 
the differential equation becomes 


on 2 ae ” 
Vive — ce} + nV, + (l-é Wp, 0 (4.12) 
with the conditions 
¢ =0 when n=] 
n 
(4.13) 
¢ ~% En when n> 


If (4.11) is substituted into (4.12) and powers of n collected, there 
results a series of differential equations of the form 


Ly = B,U, > Lo = 2BaUo > Ly = nBU, — Py? n> 3 (4.14) 


where Ly is the Legendre operator L, (1 — ¢*)F" — 2cF + n(n — lites 
The first two of these can be solved immediately, but the important point 
is that 


2 ' 
Pp," 8 I (1 = PU, Figg ee tm + 1m 1 PUL Fagg op 


(4.15) 


al 


so that Fa will depend only on lower orders. 
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Remembering that the Legendre functions of second kind Q,, (é) have 
natural unwanted singularities at € = +1 forall n, the desired homo- 
geneous solutions of (4.14) are 


anPy-1 (5) 


where the a, are arbitrary constants. Furthermore, a particular solution 


of Ly = nB OU, is ~BT hy» and it is noteworthy that both contributions are 
polynomials of relatively simple structure in €&€. Now it can be shown by 
induction 3 that the equation Ly = Pio has a particular solution Qn~2 
(say) which is a polynomial of degree n-2 when p,_, is itself a poly- 


nomial —as jis obvious for the first few terms in (4.14). In fact, a backward 


recursive algorithm to compute the coefficients of qn-2 (5) from those of 
Py (é) can be set up which supplies both the computational mechanism and 
desired proof of existence at the same time (the important consideration is 
clearly whether or not a series assumption for qn-2 (6) leads to termina- 
tion). 


On this basis, the complete sequence can be written 


(4.16) 


F, * oP ke) ~ 87 c) = 4 _o(e), n>2 


and the only remaining task is determination of the oy, constants from 
the boundary conditions. 


4.4 Boundary Conditions 

In 4.3 the general solution structure was described, but n=0 isa 
special case which does not contain any constants. It can be shown that 
Fy = — so that the leading term in the expression for » is €&n and, 
therefore, ¢~ En as n+ which means that (4.13) is satisfied. 


With regard to the surface boundary condition (4.13), the constants 


a, must be chosen. Consider for instance, the second-order approximation 


bo = nét (a — BT, I/n + Cané — BT I/né 


Sik A et a 
a A A A ral a al 
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which can also be written in terms of cosw rather than the Chebyshev 
polynomials as 


B49 Bo 
bo = a,/n + Seas cosw — —F cos2w 
n n 


giving 


I) 

Ss = -a, + (1 + B, — 2a) cosw — 2B, cos2w 
n=] 

‘| Obviously ay and do cannot be chosen to make this vanish for all w 

unless it happens that B, is zero. If B, = 0 then we choose a, = 0, 


Oy = (1 + B, )/2 to give 


B, 1+ B, 
do = [n —~— + —y— ] Cos (4.17) 
n 2n 


and if Bo #0 then the term in cos2w will not satisfy the boundary 
condition. Ellipsoids are governed by choice of B, in the range 

(-1 < B < 1) but (4.17) is not the solution for arbitrary ellipsoids 
even though the B, (n > 2) vanish for both. The only exact case for 
(4.17) is By = 0 which generates a circular profile and gives 


$5 = ( 2 cosw 


2n 


as the correct potential for a sphere. Successive approximations for true 
ellipsoids are obtained by taking more terms and, in fact, 6 has been 

| worked out by hand, but this academic topic will not be pursued further 

| here. More details can be obtained from ref. 3. 

| 


| However, the procedure described above is quite general. If the ser‘es 
for is terminated at the n" level, then the By term will not satisfy 
the surface boundary condition. For all higher-order terms, the a, are 
determined by a set of n-2 linear simultaneous equations readily solved by 
Z standard methods. Since the number of terms required for bodies of practical 
interest is far less than the number of points required to describe the 


surface shape, the solution of this system of equations is a trivial task 
by modern computing standards. Typical term and point numbers are illustrated 
by the examples given later. 


4.5 Velocity Calculation 


For engineering purposes the quantity of greatest interest is the 
surface velocity which is given by 


36, /3w 
9, * - “a57a ae 
where s_ is the arc length measured from the trailing edge. There is a 


need for some care in using this formula since it is clear from (4.4) that 


) -t/t 


7 1g (4.19) 


which vanishes at the trailing edge (w= 0). However, it is also clear that 
¢ can be expressed in the form 


@ = G(n) + G)(n) cosw + G(n) cos2w +... 


once the a, constants have been found. Therefore, 


2¢ = —[G)(n) sinw + 26,(n) sindw + ...] 


CIN) 


so that a factor sinw can be extracted to give 


ao = -sinu[G, + 26,U, + ...] = -sinwG(e,n) (4.20) 


where the U, are again Chebyshev polynomials of the second kind. Equa- 
tion (4.20) shows that the sinw factor can be used to cancel the zero in 
the denominator of (4.18) with the final, computationally sound result 


_ [2 sin(w/2)]™/™ cos(w/2) G(E,1 , 
q, = Le sinle/2ii"/* easto/2) 6(6.1) ai) 


On this basis, direct computation of velocity given either the shape itself 
or B. coefficients is quite straightforward. The central operation is 

a subroutine which calculates the coefficients of the polynomial powers of 
each succeeding Fa by recurrence relations and the pq algorithm 
mentioned in §4.4. These programs are described in more detail later, 
however, it is worth describing at this point the somewhat different 


philosophy behind the use of this theory as a design tool. 


The obvious iterative procedure starting from a given (desired) Qvs s 


as input,is to guess Bh coefficients and compute the axisymmetric part 
36,/8w together with ds/dw. One then argues that 36,/ dw is less sensitive 
to the coefficients and computes the new estimate of ds/dw by the algorithm 


(<) tee (26/2), 
dw n+] a, 


where Q. is the given velocity. It is then a standard procedure using the 
Fast Fourier Transform subroutine to calculate the a, coefficients (see 
e.g. 3). 


Unfortunately, the presence of the singular factor (2 sin(w/2)7/" 


in the velocity formula makes this computationally unacceptable except for 
cusped bodies += 0, and even then the front stagnation region w® rt 
causes difficulties. Consequently a procedure based on the wholly non- 
singular, nonzero group 


= G/|g| 


was devised and forms the basis of the programs discussed later. The under- 
lying principle is the same. Let the input desired velocity be denoted by 
q, and quantities derived from it likewise by ah etc., then the new 
algorithm is 


i) 


Iglu = a (4.22) 


This is put into effect by using the calculated arc length from a given set 

of coefficients to interpolate for x as a function of w by dividing out 

the singular factor from the given velocity. Since G, is also known from 

the n-th level of coefficients, the new estimate for |g| and hence the new 
coefficients can be calculated from (4.22). 


4.6 Solution Procedure 


The body and its associated velocity distribution are essentially des- 
cribed by NP+1 point values in the upper half plane, shown in fig. 2, and 
taken from the trailing edge as origin. However, it is NP, not NP+l, 
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Figure 2. Body Geometry. 


that is the crucial input parameter. There is a reason for this choice. 
Because of the use of a Fast Fourier Transform algorithm which uses a fold- 


ing technique, NP must contain only the factors 2, 3 or 5. 


Based on the provision of NP+] data points, whether as (x,y)-values 
(direct mode), or as (Q,s)-values (design mode), the solution procedure 
depends upon a central subroutine which implements the conformal mapping, 
described in section 4.2, from a set of Fourier coefficients. These coef- 
ficients are computed internally from the given data. 


It is important to realize with regard to the design program that a 
given velocity does not necessarily design a closed body. The program out- 
puts the body "closest" in a least squares (or rather Fourier) sense to 
having the desired velocity by compromising the overall level of maximum 
velocity whilst attempting to preserve the details. Thus, the designer 
is forced in turn to make compromises on his input levels and repeat the 
process. 


This is an inescapable feature of design methods and has nothing to do 
with the structure of this particular program or even of axisymmetric as 
opposed to two-dimensional] theory. For further discussion, see 3. 
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5.0 DESCRIPTION OF THE POTENTIAL-FLOW COMPUTER PROGRAM 


5.1 Input 


The input to the computer program consists of four types of cards. Card 
type 1 contains the title of the flow problem under consideration. 


Lala] «fs felafelohghstdneh hehe cesta pees) bibahbadadelebopel ip 474 4 oy rpep rho 
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TITLE 


Load Sheet for Card 1 


The title is input as an alphanumeric array with a format of (18A4). 


Card type 2 requires the following control data to be specified. The 
input format is (715). 


HRSABARCEDEIoS TSI Re SSeS Seeeee gens 
T p H 
pn ef wtf ae Pte 


f Load Sheet for Card 2 


f NA Number of A, coefficients to be used (for Fourier transforms ) 
i (default is NP). 
F 
i NP Number of surface elements on half of the body from trailing edge 
i to leading edge. 
F | 
NT Number of terms to be used in the axisymmetric series (default is 21). 
E NTP 0 or 1 (formerly used for trailing edge type specification, but 
i currently not operative). 
i NC Number of iterative cycles to be used in the axisymmetric inverse 
i calculation (default is 20). 
p 

IQPT Option for selective runs. 


direct axisymmetric velocity, then uses the computed velocity 
and normalized arc lengths to calculate the inverse shape 


IMPT=0 program uses the body shape input (x,y data), calculates the 
| (both direct and inverse calculations). 


I9PT=1 program calculates only the direct flow solution from the 
input x,y data. 

I9@PT=2 program reads in the user's input velocity vs normalized 
arc-length data and computes the inverse shape only. 


IPUNCH Option for punch output. 
IPUNCH=1 program will generate the x,y data and the Q,S (velocity 
and normalized arc length) data. Two points per card. 
IPUNCH=0 no punch output is desired. 


Card type 3 contains the following body shape information to be specified. — 
The input is in (2F20.10) format. (This input is required only if I@PT=0 or 1.) 


fafa] «fsfel fe] poh sehsh fea papa posh passbep subpart 


Load Sheet for Card 3 


X(J) x-coordinate of the input axisymmetric body. 
Y(J) y-coordinate of the input axisymmetric body. 
Note: The index J should start from 1 to (NP+1) where NP is 


explained in card type No. 2. One pair of x,y per card up to 
(NP+]1) pairs. The order of input is from trailing edge to leading 
edge. The trailing-edge point should always be 0, so that the x's 
go from 0 to negative value at leading edge. Only the upper half 
of the body is needed to be input as illustrated in Figure 2. 
Card type 4 contains the velocity data to be specfiied. (This input is 


required only if IMPT=2.) The format for this input is 2F20.10. 


TL Tobe badd shell bebsh-bdude pop pspspepeprhe 
QA(J) SA(J) 


Load Sheet for Card 4 


QA(J) velocity at point J. 
SA(J) normalized arc length at point J. 
Note: One pair of QA vs SA per card. Subscript J starts at 1] and goes 


to (NP+1). There should be (NP+1) pair of (QA,SA) cards to be input. 
The order of input is from trailing edge to leading edge. 
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5.2 Qutput 


The output from the inviscid potential-flow program first consists of 
a repeat of the input data used for the calculation. These are NA, the 
number of Fourier coefficients used; NTP, a parameter previously used 
to specify the trailing-edge type, but no longer used in the present code; 
NP, the number of body elements; NC, the number of iteration cycles 
specified (valid for inverse calculations only); NT, the number of terms 
in the axisymmetric expansion; IMPT, the parameter indicating direct or 
inverse mode. The specified body coordinates (direct case) or velocity and 
arc length (inverse case) are also given in tabular form. 


For the direct calculation, the calculated output is presented in 
tabular form. The columnar headings are J, OM(J), A(J), B(J), S(J), P(J), 
Q(J). These are: 


J designates the body point (last J=NP+1) 


QM(J) angular displacement of body points measured from trailing 
edge (0 < OM(J) < 7m) 


A(J) x-location of body point 

B(J) y-location of body point 

S(J) normalized arc length of body point 

P(J) an internally used special axisymmetric derivative 

Q(J) computed value of the velocity 

For the inverse calculation a running count of the iteration cycle and 
computed error is printed. If this value of the error falls below 107° the 
iterations are stopped and final results are printed. If this criterion is 


not satisfied, the iteration proceeds for NC cycles and the values after 
NC iterations are printed. 


The final tabular results have the same headings in the inverse case 
as in the direct case, i.e., J, OM(J), A(J), Bld), S(J), P(J), Q(J), but 
only J, OM(J) and P(J), have the same meaning. The others are 


A(J) calculated x-location of body point 
B(J) calculated y-location of body point 


S(J) normalized arc length of calculated body 


Q(J) "least squares" velocity fit at body points, for converged 
shape (c.f. section 4.6). 


5.3 Sample Calculations 


Altogether, five bodies were used in the test cases described below. 
The first three, Bl, B2 and B3, were characterized by having known Fourier 
coefficients for the mapping. Two more bodies, K1 and K2, were used, and 
these were chosen from boundary-layer design considerations. Therefore, as 
far as the theory described here is concerned, they were quite arbitrary, 
being defined only by (x,y) points. Figure 3 shows all five bodies. 


Figure 4 shows the results of a direct calculation on bodies Bl, B2, 
B3. The velocities are compared with the higher-order Neumann program! 9 ) 
and it can be seen that the numbers of terms used in the present method were 
adequate to give excellent results. 


In case there was a hidden bias toward the first three bodies just 
because they were originally derived from a conformal mapping, the direct 
solution procedure was applied to bodies KI and K2. Figure 5 shows, once 
again, excellent results for the velocity distribution when compared with 
higher-order Neumann. Note that because these were thought to be harder 
cases (e.g., Kl has a flat portion as can be seen from Figure 5, and this 
is a poorly convergent feature) more terms were used (NA=21). Note also that 
Kl and K2 have nonintegral values of t/r1. 


Finally, Figure 6 shows the results of using the program in the inverse 
(design) mode. Only the harder cases, Kl and K2, were used and it was 
expected that the original shapes should be returned. As can be seen in the 
figure, they were. 


test bodies. 


The five 
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Figure 3. 
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Figure 4. Velocity distribution for bodies Bl, B2, 83 using oresent method 
and higher-order Neumann 
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Figure 4. 


HIGHER-ORDER NEUMANN 
* = FOURIER COEFFICIENTS COMPUTED 


Figure 5. Velocity distribution for bodies Kl, K2 using present method 
and higher-order Neumann. 
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Figure 6. Computed bodv shapes for Kl and K2 using the (0,S) inputs 
from figure 5. 
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